library(tidyverse)
library(readxl)
library(mgcv)
library(emmeans)
library(ggeffects)
library(lme4)
library(DT)
library(ggpubr)
library(cowplot)
source("scripts/utils/utils.R")
modernacolors()## <environment: R_GlobalEnv>
# load ELISpot data
temp <- read_excel("processed_data/ELISpot.xlsx")
dp <- temp[!is.na(temp$counts), ] %>%
filter(`low cell number (<0.15M)`== "no",
day != "143")%>%
mutate(counts = counts + 0.001)
d <-
dp %>%
filter(interval != "PBS")%>%
mutate(animal_id = factor(animal_id),
interval = factor(interval),
day = factor(day))dp %>%
filter(interval != "PBS") %>%
mutate(interval = factor(interval, levels = c("PBS", 0:8))) %>%
mutate(
interval = forcats::fct_recode(
interval,
"8-week interval" = "8",
"6-week interval" = "6",
"4-week interval" = "4",
"3-week interval" = "3",
"2-week interval" = "2",
"1-week interval" = "1",
"Prime only" = "0"
)
) %>%
mutate(day = forcats::fct_recode(factor(day), "day 64" = "64")) %>%
ggplot(aes(x = interval, y = counts)) +
geom_boxplot(outlier.shape = NA) +
scale_shape_manual(values = c(16, 5)) +
scale_y_log10(breaks = c(100, 500, 1000, 1500)) +
geom_jitter(aes(x = interval, y = counts)) +
stat_summary(
fun.min = mean,
fun = mean,
fun.max = mean,
color = "red",
width = .5,
geom = "errorbar",
size = 1
) +
facet_wrap(tissue ~ coating, scales = "free") +
theme_bw() + theme(panel.grid.major = element_line(colour = "grey100")) +
labs(
x = "",
y = "antigen specific \n ASC counts/million cells",
color = "",
linetype = "",
title = "Antigen specific spleen and BM results"
) +
theme(
title = element_text(size = 14),
axis.text.x = element_text(
angle = 45,
hjust = 1,
vjust = 1
),
axis.title = element_text(size = 14),
strip.text = element_text(size = 14),
legend.text = element_text(size = 14),
legend.title = element_text(size = 14),
plot.title = element_text(size = 14, face = "bold"),
legend.position = "bottom"
)#-------------------------------------------------------------------------------
# Fit the model for spleen data
#-------------------------------------------------------------------------------
# boxcox function to fit boxcox regression
boxcox_fit <- function(dat) {
require(MASS)
if (length(unique(dat$day)) == 1) {
out <- boxcox(counts ~ interval, data = dat)
optimal <- out$x[which.max(out$y)]
bctran <- make.tran("boxcox", optimal)
lm_fit <- with(bctran,
lm(linkfun(counts) ~ interval, data = dat))
} else{
out <- boxcox(counts ~ interval * day, data = dat)
optimal <- out$x[which.max(out$y)]
bctran <- make.tran("boxcox", optimal)
lm_fit <- with(bctran,
lm(linkfun(counts) ~ interval * day, data = dat))
}
return(list(lm_fit = lm_fit, lambda = optimal))
}
# box cox model fit for spleen data
d_tissue <- subset(d, tissue == "spleen") %>%
filter(!(coating == "NTD" & interval == "0"))
boxcox_fits <-
d_tissue %>%
split(.$coating) %>%
map( ~ boxcox_fit(.x))lm_fits <-
purrr::transpose(boxcox_fits)$lm_fit
# obtain optimal lambda for boxcox transformation
lambdas <-
purrr::transpose(boxcox_fits)$lambda %>% unlist()
assumption_plots <-
imap(
lm_fits ,
~ assumption_check(.x) %>%
plot_grid(
ggdraw() + draw_label(
paste0("Spleen ", .y, "-specific model residual diagnostics "),
x = 0.03,
hjust = 0
),
.,
ncol = 1,
rel_heights = c(0.1, 1),
align = "v"
)
)
# obtain emmeans
lm_grid <- ref_grid(lm_fits$S2P, cov.keep = c("interval", "day"))
emmt_dat <-
emmeans(lm_grid, specs = ~ interval, type = "response") %>%
as_tibble()ggarrange(plotlist = assumption_plots)#-------------------------------------------------------------------------------
# Figure 4A. Spike-specific Antibody Secreting Cells and Long-Lived Plasma Cells.
#-------------------------------------------------------------------------------
# fold change figure
foldchange1 <-
emmt_dat %>%
mutate(
interval = forcats::fct_recode(
interval,
"8-week interval" = "8",
"6-week interval" = "6",
"4-week interval" = "4",
"3-week interval" = "3",
"2-week interval" = "2",
"1-week interval" = "1",
"Prime only" = "0"
)
) %>%
mutate(facet = "1 week post-boost") %>%
ggplot(aes(x = interval , y = response, color = interval)) +
geom_point(position = position_dodge(width = 0.5), size = 2) +
geom_errorbar(
aes(ymin = lower.CL, ymax = upper.CL),
width = .1,
position = position_dodge(width = 0.5),
size = 1
) +
scale_color_manual(
values = c(
modernadarkgray,
modernadarkblue2,
modernaorange,
modernapurple,
modernagreen,
modernablue,
modernared
)
) +
theme_bw() + theme(panel.grid.major = element_line(colour = "grey100")) +
scale_y_log10(breaks = c(100, 500, 1000, 1600)) +
#ggtitle("S2P specific ASCs") +
facet_grid( ~ facet) +
labs(x = "", y = "Estimated counts/million cells") +
geom_hline(yintercept = emmt_dat$lower.CL[c(6, 7)],
color = "darkcyan",
lty = "dashed") +
theme(
axis.text = element_text(size = 12),
axis.text.x = element_text(
angle = 45,
hjust = 1,
vjust = 1
),
axis.title = element_text(size = 14),
strip.text = element_text(size = 14),
legend.text = element_text(size = 14),
legend.title = element_text(size = 14),
legend.position = "None",
plot.title = element_text(size = 14, face = "bold"),
panel.spacing = unit(0, "lines")
)
# boxplot
box1 <-
dp %>%
filter(interval != "PBS", tissue == "spleen", coating == "S2P") %>%
mutate(interval = factor(interval, levels = c("PBS", 0:8))) %>%
mutate(facet = "1 week post-boost") %>%
mutate(
interval = forcats::fct_recode(
interval,
"8-week interval" = "8",
"6-week interval" = "6",
"4-week interval" = "4",
"3-week interval" = "3",
"2-week interval" = "2",
"1-week interval" = "1",
"Prime only" = "0"
)
) %>%
mutate(day = forcats::fct_recode(factor(day), "day 64" = "64")) %>%
ggplot(aes(x = interval, y = counts, color = interval)) +
geom_boxplot(outlier.shape = NA) +
scale_shape_manual(values = c(16, 5)) +
scale_y_log10(breaks = c(100, 500, 1000, 1500)) +
geom_jitter(aes(x = interval, y = counts)) +
stat_summary(
fun.min = mean,
fun = mean,
fun.max = mean,
color = "grey",
width = .5,
geom = "errorbar",
size = 1
) +
facet_grid( ~ facet) +
theme_bw() + theme(panel.grid.major = element_line(colour = "grey100")) +
labs(
x = "",
y = "Counts/million cells",
color = "",
linetype = "",
title = "S2P specific ASCs"
) +
scale_color_manual(
values = c(
modernadarkgray,
modernadarkblue2,
modernaorange,
modernapurple,
modernagreen,
modernablue,
modernared
)
) +
theme(
axis.text = element_text(size = 12),
title = element_text(size = 14),
axis.text.x = element_text(
angle = 45,
hjust = 1,
vjust = 1
),
axis.title = element_text(size = 14),
strip.text = element_text(size = 14),
legend.text = element_text(size = 14),
legend.title = element_text(size = 14),
plot.title = element_text(size = 14, face = "bold"),
legend.position = "None",
)
cowplot::plot_grid(
plotlist = list(box1, foldchange1),
nrow = 2,
ncol = 1,
labels = c("A", "B")
)#-------------------------------------------------------------------------------
# Supplementary Table 3
#-------------------------------------------------------------------------------
emm_elispot <-
emmeans(ref_grid(lm_fits$S2P, cov.keep = c("interval", "day")),
specs = ~ interval,
type = "response")
set.seed(100)
sum_elispot <-
contrast(
regrid(emm_elispot, transform = "log10"),
"pairwise",
infer = TRUE,
adjust = "mvt",
ratio = FALSE,
level = .95
) %>% as_tibble() %>%
bind_rows()
sum_elispot %>%
as_tibble() %>%
dplyr::select(contrast, estimate, p.value) %>%
filter(p.value < 0.05) %>%
mutate(estimate = 10 ^ estimate) %>%
mutate(contrast = str_replace(contrast, " - ", " / ")) %>%
rename(
"Pairwise comparison" = contrast,
"Fold change" = "estimate",
"Adjusted P-value" = "p.value"
) %>%
DT::datatable(caption = "Statistical comparisons of S2-P–specific ASCs 1 week following dose 2 between mRNA-1273 (10 μg) dosing interval") %>%
formatSignif(columns = c('Fold change', 'Adjusted P-value'),
digits = 7)# Supplementary Table 3 visualization
t3v <-
sum_elispot %>%
as_tibble() %>%
separate(contrast, c("a", "b"), sep = " - ") %>%
mutate(a = str_remove_all(a , "interval"),
b = str_remove_all(b , "interval")) %>%
mutate(
a = case_when(
str_detect(a, "1") ~ paste(a, "week"),
str_detect(a, "0") ~ str_replace(a, "0", "Prime only"),
TRUE ~ paste(a, "weeks")
),
b = case_when(
str_detect(b, "1") ~ paste(b, "week"),
str_detect(b, "0") ~ str_replace(b, "0", "Prime only"),
TRUE ~ paste(b, "weeks")
)
) %>%
mutate(contrast = paste(a, "/", b)) %>%
mutate(day = "1 week post-boost)") %>%
rename(
"Pairwise comparison" = contrast,
"Fold change" = "estimate",
"Adjusted P-value" = "p.value",
) %>%
mutate(`Pairwise comparison` = str_remove_all(`Pairwise comparison` , "interval")) %>%
ggplot(aes(x = `Pairwise comparison` , y = `Fold change`)) +
geom_point(size = 3) +
geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
width = .1,
#position = position_dodge(width = 2),
size = 1) +
scale_color_manual(
values = c(
"#ABCFE5",
"#93C4DE",
"#78B5D8",
"#60A6D1",
"#4A97C9",
"#3787C0",
"#2575B7",
"#1664AB",
"#09539D",
"#084185"
)
) +
theme_bw() + theme(panel.grid.major = element_line(colour = "grey100")) +
scale_y_continuous(breaks = log10(c(0.1, 0.2, 0.5, 1, 2, 3)),
labels = c("1/10", "1/5", "1/2", 1, 2, 3)) +
facet_grid( ~ day) +
geom_hline(yintercept = 0,
color = "darkcyan",
lty = "dashed") +
labs(title = "S2-P-specific ASCs",
y = "Estimated Fold-Change",
x = "Prime-boost dosing interval") +
coord_flip() +
theme(
axis.text = element_text(size = 12),
axis.title = element_text(size = 12),
strip.text = element_text(size = 10),
legend.text = element_text(size = 14),
legend.title = element_text(size = 14),
plot.title = element_text(size = 14, face = "bold"),
panel.spacing = unit(0, "lines")
)#-------------------------------------------------------------------------------
# Fit the model for BM data
#-------------------------------------------------------------------------------
# BM
d_bm <- subset(d, tissue == "BM")
boxcox_fits <-
d_bm %>%
split(.$coating) %>%
map(~ boxcox_fit(.x))lm_fits <-
purrr::transpose(boxcox_fits)$lm_fit
# obtain emmeans
lm_grid <- ref_grid(lm_fits$S2P, cov.keep = c("interval", "day"))
emmt_dat <-
emmeans(lm_grid, specs = ~ interval | day, type = "response") %>%
as_tibble()assumption_plots <-
imap(
lm_fits ,
~ assumption_check(.x) %>%
plot_grid(
ggdraw() + draw_label(
paste0("BM ", .y, "-specific model residual diagnostics "),
x = 0.03,
hjust = 0
),
.,
ncol = 1,
rel_heights = c(0.1, 1),
align = "v"
)
)
ggarrange(plotlist = assumption_plots)#-------------------------------------------------------------------------------
# Figure 4B. Spike-specific Antibody Secreting Cells and Long-Lived Plasma Cells.
#-------------------------------------------------------------------------------
foldchange2 <-
emmt_dat %>%
mutate(
interval = forcats::fct_recode(
interval,
"8-week interval" = "8",
"6-week interval" = "6",
"4-week interval" = "4",
"3-week interval" = "3",
"2-week interval" = "2",
"1-week interval" = "1",
"Prime only" = "0"
)
) %>%
mutate(day = forcats::fct_recode(
factor(day),
"4 week post-boost" = "87",
"24 week post-boost" = "226"
)) %>%
mutate(day = factor(day, c("4 week post-boost" , "24 week post-boost"))) %>%
ggplot(aes(x = interval , y = response, color = interval)) +
geom_point(position = position_dodge(width = 0.5), size = 2) +
geom_errorbar(
aes(ymin = lower.CL, ymax = upper.CL),
width = .1,
position = position_dodge(width = 0.5),
size = 1
) +
theme_bw() + theme(panel.grid.major = element_line(colour = "grey100")) +
scale_y_log10(breaks = c(100, 500, 1000, 1600)) +
labs(x = "", y = "Estimated counts/million cells") +
geom_hline(yintercept = emmt_dat$lower.CL[c(6)],
color = "darkcyan",
lty = "dashed") +
geom_hline(yintercept = emmt_dat$lower.CL[c(13)],
color = "darkcyan",
lty = "dashed") +
facet_grid(~ day) +
scale_color_manual(
values = c(
modernadarkgray,
modernadarkblue2,
modernaorange,
modernapurple,
modernagreen,
modernablue,
modernared
)
) +
theme(
axis.text = element_text(size = 12),
axis.text.x = element_text(
angle = 45,
hjust = 1,
vjust = 1
),
axis.title = element_text(size = 14),
strip.text = element_text(size = 14),
legend.text = element_text(size = 14),
legend.title = element_blank(),
legend.position = "None",
plot.title = element_text(size = 14, face = "bold"),
panel.spacing = unit(0, "lines")
)
box2 <-
dp %>%
filter(interval != "PBS", tissue == "BM", coating == "S2P") %>%
mutate(interval = factor(interval, levels = c("PBS", 0:8))) %>%
mutate(
interval = forcats::fct_recode(
interval,
"8-week interval" = "8",
"6-week interval" = "6",
"4-week interval" = "4",
"3-week interval" = "3",
"2-week interval" = "2",
"1-week interval" = "1",
"Prime only" = "0"
)
) %>%
mutate(day = forcats::fct_recode(
factor(day),
"4 week post-boost" = "87",
"24 week post-boost" = "226"
)) %>%
mutate(day = factor(day,
levels = c("4 week post-boost", "24 week post-boost"))) %>%
ggplot(aes(x = interval, y = counts, color = interval)) +
geom_boxplot(outlier.shape = NA) +
scale_shape_manual(values = c(16, 5)) +
scale_y_log10(breaks = c(100, 500, 1000, 1600)) +
geom_jitter(aes(x = interval, y = counts)) +
stat_summary(
fun.min = mean,
fun = mean,
fun.max = mean,
color = "grey",
width = .5,
geom = "errorbar",
size = 1
) +
scale_color_manual(
values = c(
modernadarkgray,
modernadarkblue2,
modernaorange,
modernapurple,
modernagreen,
modernablue,
modernared
)
) +
theme_bw() + theme(panel.grid.major = element_line(colour = "grey100")) +
labs(
x = "",
y = "Counts/million cells",
color = "",
linetype = "",
title = "S2-P-specific LLPCs"
) +
facet_grid(~ day) +
theme(
axis.text = element_text(size = 12),
title = element_text(size = 14),
axis.text.x = element_text(
angle = 45,
hjust = 1,
vjust = 1
),
axis.title = element_text(size = 14),
strip.text = element_text(size = 14),
legend.text = element_text(size = 14),
legend.title = element_text(size = 14),
legend.position = "None",
plot.title = element_text(size = 14, face = "bold")
)
cowplot::plot_grid(
box1,
box2,
foldchange1,
foldchange2,
labels = c("A", "C", "B", "D"),
align = 'hv',
rel_widths = c(1, 2),
rel_heights = c(1, 1)
)#ggsave(filename = "figures_tables/ELISpot/ELISpot_boxplot_F4.png", width = 14, height=10)
#ggsave(filename = "figures_tables/ELISpot/ELISpot_boxplot_F4.pdf", width = 14, height=10)#-------------------------------------------------------------------------------
# Supplementary Table 4
#-------------------------------------------------------------------------------
emm_elispot_BM <-
emmeans(ref_grid(lm_fits$S2P, cov.keep = c("interval", "day")),
specs = ~ interval |
day,
type = "response")
set.seed(100)
emm_elispot_BM <-
contrast(
regrid(emm_elispot_BM , transform = "log10"),
"pairwise",
infer = TRUE,
adjust = "mvt",
ratio = FALSE,
level = .95
) %>% as_tibble() %>%
bind_rows()
emm_elispot_BM %>%
as_tibble() %>%
dplyr::select(day, contrast, estimate, p.value) %>%
filter(p.value < 0.05) %>%
mutate(day = as.numeric(as.character(day))) %>%
arrange(day) %>%
mutate(estimate = 10 ^ estimate) %>%
mutate(contrast = str_replace(contrast, " - ", " / ")) %>%
rename(
"Pairwise comparison" = contrast,
"Fold change" = "estimate",
"Adjusted P-value" = "p.value",
"Day" = "day"
) %>%
DT::datatable(
caption = "Statistical comparisons of S2-P–specific LLPCs 4 weeks following
dose 2 (Day 87) and 24 weeks following dose 2 (Day 227) between mRNA-1273
(10 μg) dosing intervals"
) %>%
formatSignif(columns = c('Fold change', 'Adjusted P-value'),
digits = 7)t4v <-
emm_elispot_BM %>%
as_tibble() %>%
separate(contrast, c("a", "b"), sep = " - ") %>%
mutate(a = str_remove_all(a , "interval"),
b = str_remove_all(b , "interval")) %>%
mutate(
a = case_when(
str_detect(a, "1") ~ paste(a, "week"),
str_detect(a, "0") ~ str_replace(a, "0", "Prime only"),
TRUE ~ paste(a, "weeks")
),
b = case_when(
str_detect(b, "1") ~ paste(b, "week"),
str_detect(b, "0") ~ str_replace(b, "0", "Prime only"),
TRUE ~ paste(b, "weeks")
)
) %>%
mutate(contrast = paste(a, "/", b)) %>%
mutate(day = factor(day, levels = c(87, 226))) %>%
mutate(day = forcats::fct_recode(
day,
"4 week post-boost" = "87",
"24 week post-boost" = "226",
)) %>%
rename(
"Pairwise comparison" = contrast,
"Fold change" = "estimate",
"Adjusted P-value" = "p.value",
) %>%
mutate(`Pairwise comparison` = str_remove_all(`Pairwise comparison` , "interval")) %>%
ggplot(aes(x = `Pairwise comparison` , y = `Fold change`)) +
geom_point(size = 3) +
geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
width = .1,
#position = position_dodge(width = 2),
size = 1) +
scale_color_manual(
values = c(
"#ABCFE5",
"#93C4DE",
"#78B5D8",
"#60A6D1",
"#4A97C9",
"#3787C0",
"#2575B7",
"#1664AB",
"#09539D",
"#084185"
)
) +
theme_bw() + theme(panel.grid.major = element_line(colour = "grey100")) +
scale_y_continuous(breaks = log10(c(0.1, 0.2, 0.5, 1, 2, 3)),
labels = c("1/10", "1/5", "1/2", 1, 2, 3)) +
geom_hline(yintercept = 0,
color = "darkcyan",
lty = "dashed") +
facet_grid(~ day) +
labs(title = "S2-P-specific LLPC",
# subtitle = "4 weeks following dose 2 (Day 87) and 24 weeks following dose 2 (Day 227)
# between mRNA-1273 (10 μg) dosing intervals",
y = "Estimated Fold-Change",
x = "Prime-boost dosing interval") +
coord_flip() +
theme(
axis.text = element_text(size = 12),
axis.title = element_text(size = 12),
strip.text = element_text(size = 10),
legend.text = element_text(size = 14),
legend.title = element_text(size = 14),
plot.title = element_text(size = 14, face = "bold"),
panel.spacing = unit(0, "lines")
)
cowplot::plot_grid(
plotlist = list(t3v, t4v),
labels = c("A", "B"),
nrow = 1,
ncol = 2,
rel_widths = c(0.5, 1)
)ggsave(filename = "figures_tables/ELISpot/emmeans_BM_spleen.png",
width = 16,
height = 8)
ggsave(filename = "figures_tables/ELISpot/emmeans_BM_spleen.pdf",
width = 16,
height = 8)sessionInfo()## R version 3.6.3 (2020-02-29)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Catalina 10.15.7
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] MASS_7.3-51.5 cowplot_1.1.1 ggpubr_0.4.0 DT_0.24
## [5] lme4_1.1-28 Matrix_1.2-18 ggeffects_1.1.3 emmeans_1.7.5
## [9] mgcv_1.8-31 nlme_3.1-144 readxl_1.4.0 forcats_0.5.1
## [13] stringr_1.4.0 dplyr_1.0.9 purrr_0.3.4 readr_2.1.2
## [17] tidyr_1.2.0 tibble_3.1.8 ggplot2_3.3.6 tidyverse_1.3.2
##
## loaded via a namespace (and not attached):
## [1] fs_1.5.2 lubridate_1.8.0 httr_1.4.3
## [4] tools_3.6.3 backports_1.4.1 bslib_0.4.0
## [7] utf8_1.2.2 R6_2.5.1 DBI_1.1.3
## [10] colorspace_2.0-3 withr_2.5.0 tidyselect_1.1.2
## [13] compiler_3.6.3 textshaping_0.3.6 cli_3.3.0
## [16] rvest_1.0.2 xml2_1.3.3 labeling_0.4.2
## [19] sass_0.4.2 scales_1.2.0 mvtnorm_1.1-3
## [22] systemfonts_1.0.2 digest_0.6.29 minqa_1.2.4
## [25] rmarkdown_2.14 pkgconfig_2.0.3 htmltools_0.5.3
## [28] highr_0.9 dbplyr_2.2.1 fastmap_1.1.0
## [31] htmlwidgets_1.5.4 rlang_1.0.4 rstudioapi_0.13
## [34] jquerylib_0.1.4 generics_0.1.3 farver_2.1.1
## [37] jsonlite_1.8.0 crosstalk_1.2.0 car_3.1-0
## [40] googlesheets4_1.0.1 magrittr_2.0.3 Rcpp_1.0.9
## [43] munsell_0.5.0 fansi_1.0.3 abind_1.4-5
## [46] lifecycle_1.0.1 stringi_1.7.8 yaml_2.3.5
## [49] carData_3.0-5 grid_3.6.3 crayon_1.5.1
## [52] lattice_0.20-38 haven_2.5.0 splines_3.6.3
## [55] hms_1.1.1 knitr_1.39 pillar_1.8.0
## [58] boot_1.3-24 estimability_1.4.1 ggsignif_0.6.3
## [61] reprex_2.0.1 glue_1.6.2 evaluate_0.16
## [64] modelr_0.1.8 vctrs_0.4.1 nloptr_1.2.2.3
## [67] tzdb_0.3.0 cellranger_1.1.0 gtable_0.3.0
## [70] assertthat_0.2.1 cachem_1.0.6 xfun_0.32
## [73] xtable_1.8-4 broom_1.0.0 coda_0.19-4
## [76] rstatix_0.7.0 ragg_1.1.3 googledrive_2.0.0
## [79] gargle_1.2.0 ellipsis_0.3.2